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Abstract 

The paper begins with a short overview of the recent work done in the field of discontinuous reinforced 
composites, focusing on the different parameters which influence the material behavior of discontinuous 
reinforced composites, as well as the various analysis approaches undertaken. Based on this overview 
it became evident, that in order to investigate the enumerated effects in an efficient and comprehensive 
manner, an alternative approach to the computational ly intensive finite-element based micro-mechanics 
approach is required. Therefore an investigation is conducted to demonstrate the utility of utilizing the 
generalized method of cells (GMC). a semi-anahtical micromechanics-based approach, to simulate the 
elastic and elastoplastic material behavior of aligned short fiber composites. The results are compared 
with 1) simulations using other micromechanical based mean field models and finite element (FE) unit 
cell models found in the literature given elastic material behavior, as well as 2) finite element unit cell and 
a new semi-anah*tical elastoplastic shear lag model in the inelastic range. GMC is shown to definitely 
have a window of applicability when simulating discontinuous!} reinforced composite material behavior. 


1 Introduction 

Metal matrix composites (MMGs) have recehed considerable attention o\er the past 30 years due to 
their attracthe specific strength, stiffness, fatigue and thermal properties[Cl}Tie and Withers (1993)], [Kell} 
and Zweben (Eds) (2000)], [Sinclair and Grogson (1997)], [Lloyd (1994)]. MMC materials can be gener- 
icallv classified by their type of reinforcement, which consists of either continuous fiber (identified herein 
as CFMMC) or discontinuous particulate (termed herein as DRX) reinforcements. While both classes of 
materials share common features at a fundamental level (e.g., matrix/ reinforcement compatibility, interfa- 
cial properties sensitivity, architectural interaction (reinforcement size, shape and packing) and pro pert} 
optimization, at a practical le\el dramatic differences can and do exist. For instance, CFMMC typically 
have exceptional!} high strength and stiffness capabilities that are inherently strongly anisotropic; whereas, 
DRX material properties are typically more isotropic in nature with lower specific strength, yet ha\e excel- 
lent stiffness capability. Clearly, failure to appreciate and consider explicitly the appropriate reinforcement 
morphology can result in inaccurate performance assessments, as elaborated upon in the next section. The 
combination of good transverse properties, relatively low cost, high workability and reasonable increases in 
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performance over monolithic metallic alloys have made discontinuous reinforced materials the most commer- 
cially attractive system to date for man} stiffness driven applications [Chme and Withers (1993)], [Kelly and 
Zweben (Eds) (2000)], [Miracle and Maruyama (2000)]. 

Although manufacturing Issues and cost are probably the greatest impediments to the wide application of 
these advanced materials, it Is also clear that accurate design and life prediction tools are greatly needed to 
facilitate the implementation of these developing materials. Although closure has not been reached regarding 
the best models for use in design and life prediction, it has become evident that if a model Is ever to serve a 
purpose beyond that of basic research, it must fulfill several primary requirements. These include a significant 
level of accuracy on both the macro and micro scales, computational efficiency, and compatibility with the 
finite element method Fulfillment of these requirements allows a model to serve the composite developer 
by enabling quick and easy variation of composite parameters for material development and optimization 
purposes. Likewise, those who attempt to utilize composites in structural design are well served if the model 
is compatible with the finite element method The generalized method of cells (GMC), originally developed 
by Aboudi (1991, 1995) , Ls a good choice for implementation into modeling tools for advanced composites, 
given the requirements described above. Particularly, whenever microfield as well as macrofield quantities 
will be needed - as Is the case when attempting to account for inelastic material behavior in MMCs. 

GMC Ls a fullv- analytical micromechanics model for multi-phased materials with arbitrary periodic 
microstructures. It provides pseudo closed-form multiaxial constitutive equations for such materials, and 
allows straightforward implementation of physically-based viscoplastic deformation models, as well as arbi- 
trary failure and damage models for each phase. Further, recent independent advances have simplified the 
implementation of GMC as an elemental constituent material model in finite element analysis [Arnold et al, 
(1999)], and significantly increased the models computational efficiency [Pindera and Bednarcyk (1999)]. 
GMC has been implemented in the NASA Glenn Research Center comprehensive micromechanics analysis 
code, MAC/GMC [Arnold et al, (1999)]. The code has man} features that render it useful for design, de- 
formation modeling, and life prediction for a wide range of materials. These features include the ability to 
simulate general thermomechanical loading on composites whose geometries are represented by a library of 
continuous and discontinuous repeating unit cells, a library of nonisothermal elastic/ viscoplastic constitutive 
models, fatigue damage analysis, yield surface analysis, laminate analysis, and interface modeling. 

In the past, the various capabilities described above have been exercised extensively, however, the pre- 
dominant material system examined has been that of the CFMMC [Aboudi(1996)], [Arnold, et. al. (1996)], 
[Goldberg and Arnold (2000)], [Iyeretal. (2000)], [Bednarcyk and Arnold (2000a)], [LLssenden et al. (2000)], 
[Bednarcy k and Arnold (2000b)]. Consequently, the primary' objective of this current stud}' is to assess the 
applicability of GMC in predicting the elastic and inelastic (plastic) behavior of particulate reinforced com- 
posites (DRXs). The literature indicates that DRX behavior Ls significantl}' influenced by various charac- 
teristic parameters like fiber arrangements (staggered, non-staggered, cluster, random, packing type), fiber 
vxdume fraction, fiber aspect ratios, fiber-matrix stiffness ratios and reinforcement type (cylinder, particle, 
etc. ). A brief review of this Ls provided in the background section. Even though GMC is full}' capa- 
ble of modeling arbitrary multiphased composite architectures, to facilitate comparison with previous work 
reported in the literature and limit the scope of this investigation, only one type of DRX Ls considered, 
namely, the simulation of aligned short fiber composites with non-staggered cylindrical reinforcements and a 
hexagonal packing arrangement. The need for such an assessment stems from the fact that although GMC 
consists of analytical equations that are able to deal with both the homogenization and localization process 
for composite materials in a very numerically efficient way, the method’s linear displacement field assumption 
in combination with the imposition of displacement and traction boundary conditions in an average sense, 
gives rise to an inherent lack of normal and shear component coupling. This lack of coupling is such that 
application of a globed normed (or shear) stress field will produce only a local normal (or shear) stress field. 
Consequently, for cases (e.g., DRX materials) in which the load transfer mechanism between phases (e.g,, 
matrix and fibers) is shear dominated, the lack of normal-shear coupling may become quite problematic. 
In particular, in the case of an elastoplastic matrix regime placed between two short fibers, GMC’s lack of 
shear coupling requires both the fiber and matrix to cany' the same normal stress component. This tends 
to produce an overh compliant (soft) composite response. The results presented within demonstrate this 
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fact, and more importantly characterize the window of applicability for the GMC method within the context 
of discontinue ilsIy reinforced metallic composites. The introduction of an “artificial elasto plastic matrix 
material between the particulate reinforcements, Is proposed in an attempt to extend the applicability of 
GMC. Although, providing more reasonable macro behavior, such a strategy still results in incorrect mi- 
cromechanical stress fields within the unit oells, thereby limiting the usefulness of the current GMC model 
when calculating failure or fatigue analyses. 

With section 2, the paper begins by reviewing recent work done in the field of short fiber composites so 
as to re\eal the key parameters influencing DRX material behavior. Section 3 then articulates the various 
analytical and numerical approaches employed and compared in this paper. Finally, section 4 presents a 
comparison of GMC results with other analytical and numerical results within the elastic and elastoplastic 
range, respectively. 


2 Background: Factors Influencing DRX Materials 

The earh' work on the principles of strengthening in composites reinforced with discontinuous fibers, 
w hiskers or particulates carried out in the 1950 s and 1960 s is characterized by the use of relatively simple 
models that ck> not explicitly take into account the actual inclusion distribution through interaction between 
adjacent inclusions. Man} models have been employed to stud}' the effect of inclusion shape on the elastic 
moduli, as well as subsequent inelastic response, and also the mechanism of stress transfer from the matrix 
to a finite-length fiber. For instance, Eshelby s analytical solution to the problem of an arbitrarily shaped 
ellipsoidal inclusion embedded in an homogeneous material deformed by uniform tractions or displacements 
at infinity has formed the basis for calculating the effect he response of macroscopically homogeneous two- 
phase composites using a number of approaches [Eshelby ( 1957)], These approaches include the self-consistent 
scheme which neglects the inclusion/matrix interaction in calculating stress fields in the inclusion phase [Hill 
(1965)], and the Mori-Tanaka method (MTM) which takes this interaction into account in an approximate 
fashion [Mori and Tanaka (1973)]. The problem of an array of ellipsoidal inclusions with different aspect 
ratios embedded in an elastoplastic matrix can be treated using the Mori-Tanaka approach, as was done 
by [Brown and Clarke (1975)] investigating the effect of inclusion shape on work hardening of metal matrix 
composites. However, this method typically underestimates yielding and subsequent hardening effects due to 
the use of mean stress said strain fields in the matrix phase, and treats all inclusion distributions on the same 
footing so long as macroscopic homogeneity is preserved. Along similar lines, the so-called shear-lag anal} sis 
has been employed to stud} the effectiveness of short fibers as reinforcement using the strength-of-materiaLs 
approach to analyze stress fields around and within a finite-length fiber embedded in a surrounding matrix 
[Dow' (1963)]. While this type of analysis helps to identify shearing of the matrix as the primary mechanism 
of force transfer from the matrix to the fiber, and thus the critical fiber length over w hich the axial stress is 
introduced into the fiber from both ends, it is based on a very simplified one-dimensional anal} sis of stress 
fields w hich neglects the influence of morphology of surrounding fibers, among other things. 

Despite the relatively long history of modeling the response of discontinuous fiber composites, only 
recently have systematic in\estigations of the effects of reinforcement shape and arrangement been initiated 
for these typos of oompx>sites. Inclusion of the third dimension in the analysis of the elastic and inelastic 
response of discontinuous fiber composites increases the number of variables several fold relative to the 
analysis of continuous fiber composites. Thus, in addition to the arrangement and shape of the reinforcement 
in the plane transverse to the loading direction, variables associated with the planes parallel to the loading 
direction were included. These variables include the fiber aspect ratio, fiber spacing in the vertical and 
horizontal directions, those associated with the effect of fiber clustering, and the extent of overlap between 
adjacent columns of fibers. These added complexities typically require numerical solution procedures, such as 
the finite-element anal} sis; particular!} if complicated reinforcement shapes and arrangements are involved. 

Finite-element investigations carried out in recent years have focused on separating the effects of inclusion 
shape from the effects of inclusion packing array geometry on the overall elastic, elastoplastic and creep 
behavior. Due to man} different types of reinforcement employed in dlscontinuoush -rein forced composites, 
including whiskers, p>articulates, platelets as well as finite-length fibers, extensive analyses of inclusion shap>e 
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have been conducted for this class of composites. Examples of inclusion shapes investigated include spherical 
particles, circular cylinders with different aspect ratios, truncated circular (octagonal) cylinders, double-cone 
particles, and ellipsoidal or disk-like particles (platelets). Similarly, different packing arrays studied include 
hexagonal arrangements of inclusions in the plane of loading distributed in either regular (aligned) arrays or 
off-set (staggered) arrays in the planes parallel and perpendicular to the loading direction so as to model and 
assess the influence of distribution randomness that Is typical in DRX composites. Inclusions such as disk- 
like particles arranged in packet morphologies have also been considered, as have clustered arrangements. 
These investigations have been driven in large part by the wide range of inclusion morphologies that result 
from current material processing techniques, as well as by the wide range of shapes available for the inclusion 
phase. 

In order to reduce the complexity of a three-dimensional finite-element formulation in analyzing the 
response of discontinuous fiber, w hisker or particulate composites, various idealizations of the unit cell have 
been employed by a number of researchers (e.g. clusters of whiskers or cylindrical particles in the plane 
parallel to the applied load have been modeled using a plane strain idealization of rectangular platelets). 
This effectively reduces the problem to that of a continuously reinforced composite subjected to loading 
in the plane perpendicular to the long fiber direction. The results of such analysis for rectangular cross- 
sections with different aspect ratios in the plane of loading should be easily deduced from the early analysis 
of continuously reinforced composites with rectangularly shaped fibers [Ashton et al.(1969) ]. In the case of 
a hexagonal array of inclusions in the plane transverse to the loading direction, the problem Is often reduced 
to an axisvinmetric (i.e., twe-dimensional) problem by approximating the unit cell using a circular cylinder 
with different types of lateral boundary conditions to simulate the interaction with adjacent fibers (i.e., unit 
cells). Using this model, different fiber arrangements in the plane parallel to the applied load have been 
investigated, with varying amounts of overlap between vertical columns of adjacent fibers (e.g. [Christman 
et al. (1989)], [Tvergaard (1990)], [Dragone and Nix (I960)], [Povirk et al. (1990)], [Yang et al. (1991)], 
[Bao et al.(1991)], [Siegmundet al. (1992)]). In contrast to continuously reinforced composites, the inclusion 
arrangement in discontinuous^ reinforced composites is inherent lv three-dimensional, so that the use of 
plane models requires considerable caution. The results of full} three-dimensional finite-element analyses 
of such periodic arrays are also available (e.g. [Levy and Papazian (1990)], [Weissenbek and Rammerstorfer 
(1993)], [Abel et al.(1993)], [Weissenbek et al. (1993)]). They are, however rather restrictive in terms of 
fiber arrangements and shapes, and tend to entail considerable computational costs. Both three-dimensional 
and axisvinmetric unit cell models have been used successfully for studying the nonlinear thermomechanical 
behavior of aligned short reinforced composites. 

The literature survey presented indicates significant influence of both inclusion shape and distribution on 
the elastic, elas to plastic and creep response of discontinuous fiber composites. The extent of this influence for 
a given fiber architecture depends on the inclusion content, inclusion/matrix material property mismatch, and 
the direction of applied load with respect to the internal micro-structure. The influence of fiber architecture 
on the response of DRX composites is significant and is based on a substantial number of investigations, 
each of which w as necessarily limited in scope due to the prevalent use of the finite-element approach in 
modeling the response of the investigated architectures. Thus while a considerable body of knowledge 
has been generated that sheds light on the effectiveness of different fiber architectures in strengthening 
discontinuous fiber composites, considerably more systematic research is required to develop design guidelines 
for optimization of material performance through fiber architecture manipulation. 

However, it becomes dear, based on the above discussion, that in order to investigate the enumerated 
effects in an efficient and comprehensive manner, including parametric studies involving fiber content and 
material property mismatch variations, an alternative approach to the computation all} intensive finite- 
element based micromechanics approach is required. This is particularly true when inelastic (viscoplastic) 
theories, which typically require computationally intensive integration algorithms, are employed to model 
the response of metallic matrices. Therefore, the objective of this paper is to apply the computationally 
efficient GMC method so as to determine its suitability for predicting the macroscopic response of aligned. 
DRX materials. 
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3 Modeling Approaches 

All simulations conducted in this stud}' consist of aligned short fiber composites, w herein the particulates 
(wishers) are idealized as non-staggered. cy lindrical reinforcements, w ith a hexagonal packing arrangement. 
Similarly other basic constraints imposed are as follow s: 

• The inclusion properties are assumed isotropic and linear elastic, while the matrix is taken to be 
isotropic elastoplastic with isotropic hardening. The material data for both the inclusion and matrix 
are treated as temperature-independent and isotropic. 

• The inclusions are axisymmetric, identical in shape and size and can be characterized by an aspect 
ratio, ci r , that relate the length of an inclusion to its diameter. 

• Investigated model parameters are the fiber (inclusion) volume fraction (V/ = 0.2. 0.5). the fiber-matrix 
stiffness ratio ( s r = Ef/Em = 3. 30) and the fiber aspect ratio (a r = l/d — 5. 15, 25). 

• The inclusions and matrix are yvell bonded at their interface and remain that way during deformation. 
Thus, yve do not consider interfacial slip, fiber-matrix debonding or matrix micro-cracking. 

• No time dependence (relaxation or creep) effects are included in the present anal} sis. 

Note, obtaining experimental reference data for unidirectional short fiber composites appears to be problem- 
atic, as it has not proved passible to produce physical samples with perfectly aligned fibers. Consequently*, in 
this paper we avoid this additional complexity of orientation effect by using three-dimensional finite element 
analy sis of aligned short fiber composites as our reference standard, rather than experimental results. 


3.1 Material Data 

The material properties used for all calculations are given in Table 1. These properties are in the range 
of ty pical fiber-reinforced engineering thermo-plastics (PMCs - Material 1) and metal matrix composites 
(MMCs - Material 2), respectively. The following input data are needed for both components: Young's 
modulus E and Poisson’s ratio is. The y ield stress cry and the hardening modulus E h are needed for the 
matrix material only . 


Table 1: Prototy pical material properties used for all calculations. 



Material 1 ( 5 ^ =30) 

Material 2 (s* =3) 

Property* 

Fiber 

Matrix 

Fiber 

Matrix 

E (GPa) 

300. 

10 . 

300. 

100 . 

IS 

0.17 

0.33 

0.17 

0.33 

rry (MPa) 

- 

20 . 

- 

200 . 

Eh (GPa) 

- 

0.1 

- 

1.0 


3.2 Analytical Models 

In the past, numerous analy tical models have been examined in an attempt to determine the ‘best" 
(see e.g., [Tucker and Liang (1999)], [Aboudi (1991)]) available analytical model for predicting the elastic 
stiffness behavior of aligned short fiber reinforced composites relative to reference FEA-model represen- 
tations. Analytical models considered, were models such as: the Mori-Tanaka method (MTM)[Mori and 
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Figure 1: Sketch of inclusions-matrix shapes for the applied aligned short fiber composite models, (a) Mori- 
Tanaka model, (b) Hasliin-Sht rikman estimates of Ponte Castaneda and Willis, (c) CMC -model (a r = 5.0) 


Tanaka (1973)], self-consistent models, bounding models (e.g., \ariable constraint model (VCM) and Hashin- 
Shtrikman bounds and estimates), Halpin-Tsai equations and shear lag models (see Appendix B for details 
of the el as to plastic extension used in this stud}). For example, the recent work of [Tucker and Liang ( 1999)] 
recommended that the Mori-Tanaka model was the '‘best” choice for estimating the overall elastic stiffness 
behavior of aligned short fiber composites. 

In this stuck we will extend the scope of comparison to include the semi-analytical generalized method of 
cells model (GMC) and compare GMC with the MTM. VCM, FEM results and Hashin-Shtrikman estimates 
(HSE) 1 of [Ponte Castaneda and Willis (1995)]. Bemeniste (1987) has provided a particular!} simple and 
clear explanation of the Mori-Tanaka approach. Analytical bounds for the stiffness behavior can be obtained 
from the variable constraint model (VCM), published by [Pedersen and Withers (1992)]. From a practical 
point of view it is of interest that the lover (VCM) bound corresponds to the standard Mori-Tanaka results, 
w hile the upper (VCM) bound can also be obtained from the Mori-Tanaka approach after a so called “color 
inversion” (i.e., exchanging the roles of inclusion and matrix). These bounds obtained from the Mori-Tanaka 
approach also correspond to the Willis bounds in the case of aligned inclusions, compare for example. [Weng 
(1990)]. 

It is important to remember that the basic assumption of the Mori-Tanaka method is that the aligned 
ellipsoidal inclusions are surrounded by coaxial aligned ellipsoidal matrix domains in w hich both ellipsoids 
have the same aspect ratio (Fig. la) and that each inclusion (fiber) within a concentrated composite sees the 
average strain of the matrix. Consequently, the MTM provides only macro (average) fields and properties 
and therefore lacks the ability to address adequately problems dominated by microfield quantities, such as 
damage and inelastic behavior anal} sis. Alternatively, GMC allows both homogenization and localization of 
the associated stress and strain fields (see Appendix A), thereby providing not only macro stiflhess properties 
but also local eigenstrain field quantities. Note, that in the present investigation the GMC unit cell differs 
from MTM as different fiber-unit cell aspect ratios (Figs.2 and 3) can be realized since the w hole inclusion 
(fiber) is surrounded by a matrix material layer of constant thickness (side-to-side and end-to-end distances 
are equal) and labeled, e, subsequently. Because of this varying aspect ratio, HSE for the prediction of the 
effective moduli of the elastic inhomogeneous material given in [Ponte Castaneda and Willis (1995)] are also 
introduced (Fig. lb) to overcome the drawback of the original Mori-Tanaka method so that different fiber-unit 
cell aspect ratios can be realized Note, these Hashin-Shtrikman estimates correspond identical!} with the 
Mori-Tanaka method for the case of aligned, coaxial inclusions with equal aspect ratios of the inclusion and 
the surrounding matrix ellipsoids (see [Hu and Weng (2000)]) 

^ote. the Hashin-Shtrikman estimates differ from the Hashin-Shtrikman bounds, as [Ponte Castaneda and Willis (1995)] 
and [Willis (1977)] removed the restriction of statistically isotropic overall material symmetry (used by both HS and Willis) 
by separating the spatial distribution of inclusions from the inclusion shape. Also in the HSE. [Ponte Castaneda and Willis 
(1995)] employed only a single reference material (e.g.. the effective material like in self consistent models) instead of the two 
used by [Willis (1977)]. 
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3.3 Generalized Method of Cells 

Although, GMC Is capable of anal} zing any multiphased composite material with an arbitrary internal 
microstructure and reinforcement shape (see Appendix A), the unit cell analyzed in this stud}' is shown 
in Fig. 2. Due to the periodic boundary conditions inherent to the GMC -model the whole RYE must be 
modeled, instead of only a quarter of the unit cell (see the dashed rectangle in Fig. 2) as typically done 
in the case of FEA representations where s}mtnetric boundary conditions are applied. Apart from this 


*3 



Figure 2: Aligned short fiber GMC unit cell model, comprised of 338 subcelis. 

fact, the unit cell dimensions were chosen similar to the dimensions of the FE-model (see Fig. 3). The fiber 
was approximated with more than one subcell (sufficient for the longitudinal stiffness behavior) in order to 
reduce the error in the transverse response (i.e., the difference between E 22 and £ 33 ) and better approximate 
the targeted transversely isotropic response behavior resulting from a hexagonal packing arrangement. The 
cross sectional area of the approximated fiber is taken equal to the cross sectional area of the actual circular 
fiber. These conditions, together with the known fiber diameter, fiber volume fraction, and the distance, e 
(see Fig. 3) suffice to determine the dimensions of the unit cell for each fiber aspect ratio a r =l//d/. All 
GMC analyses were realized with the computer code MAC/GMC [Arnold et al.(1999)] developed at the 
NASA Glenn Research Center. 

3.4 Finite Element Models 

The finite element unit cell anal} zed consisted of an identical periodic, three dimensional non-staggered 
array of fibers (similar to that idealized with GMC) with the corresponding dimensions shovm in Fig. 3. The 
unit cell dimension were chosen such that b would be proportional to a (i.e., b = 3a) through the constant 
3 . In the case of hexagonal packing, 3 = \/5, and the distance between neighboring fiber ends (end-to-end 
distance, equal to c — If in Fig. 3) was set to the closest distance between two fibers (side-to-side distance, 
equal to e in Fig. 3). 

These conditions, together with the fiber diameter and volume fraction, suffice to determine the dimension 
o, b and c for the repeating vol ume element (RY E) to be analyzed. Note that a new' unit cell and its 
corresponding 3D mesh are generated for each fiber aspect ratio (l/d) and fiber volume fraction 
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GMC unit cell 


FE unit cell 




Figure 3: Example representative volume element and finite element unit cell for — 5.0 and Vp =0.2 


Stiffnesses within the elastic and plastic range of these unit cells were calculated using [ABAOUS (199S)]. 
Full}' integrated, first order, solid elements 2 were used, as this element type does not suffer from \olumetric 
locking, which is critical for accurate elastoplastic calculations. The mesh discretization containing 13.300 
elements employed is shown in Fig. 3. This mesh density is approximately three times finer than that used 
by [Weissenbek and Rammerstorfer (1993)] and was shown during the course of this work ( by convergency 
studies) to be sufficiently accurate for both elastic and inelastic analysis. Modeling of only one eighth of the 
w hole EVE (depicted on the left) was achieved due to the use of appropriate svmmetric boundary conditions. 
For example, when considering axial or transverse loading, svTnmetry requires all faces of the unit cell to 
remain plane. Consequently, to determine E\\ and v 12 the normal displacements of the front (xi = 0), left 
(X 2 = ci/2) and bottom (xj = 0) faces of the unit cell were fixed Then a master node was defined (the 
intersection point of the back, right and top face) and linear constraint equations applied on the back, right 
and top face such that the normal displacements of all nodes (except the master node) on these faces were 
set equal to the normal displacements of the master node. The tangential displacements on all faces were left 
unconstrained The average (marco) stress was then computed from the reaction force of the master node 
divided by the cross-sectional area normal to the reaction force of the investigated unit cell. Similarly, the 
average (macro) strain was computed by dividing the displacement of the master node by the corresponding 
length of the unit cell. For convenience, a one New ton load was applied to the unit cell. 

Analogous conditions were used to load the unit cell in the ^-direction to determine £22 aud v 23. 
The fifth and final independent material parameter for transversal isotropic material, the longitudinal shear 
modulus G 1 2 > could also be determined using the unit cell shown in Fig. 3, assuming the appropriate periodic 
boundary conditions were applied to each face. This more complicated analysis was not undertaken at this 
time. 


4 Results 

4.1 Elastic Material Behavior 

All of the analytical micromechanics theories used in this stud}’ predict full transversely isotropic prop- 
erties. Transverse isotropy about the x\ axis implies that the stiffness behavior is the same for any loading 

2 In [ABAQUS (1998)] this is the element known as C3D8. 
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( 1 ) 


direction in the 2-3 plane: consequently, this not only requires that £22 = £ 53 - but also that 

£22 

2(1 + ^23) 

Therefore, only five independent stiffness constants (e.g., E\\* £' 22 - /y i2* ^23 and G 12 ) need be determined for 
the case of transverse isotropy. Periodic unit cells with hexagonal packing should be transversely isotropic 
as there are six planes of symmetry. However, since periodic boundary conditions were not applied to the 
FE unit oell, we were unable to obtain G i > from the FE in\estigations. Consequently, we will limit our 
comparison for the \arious approaches to the prediction of Eu- £ 22 , "12 and is 23 only. 

Results of the normalized axial modulus (xi- direction, scaled by the matrix modulus) obtained from the 
VCM-model (MTM-model), Hashin-Shtrikman estimates, GMC-, and FE-model are shown in Fig. 4 (DRX 
like materials with a lower stiflhess ratio, e.g.. s r = 3) and Fig. 5 (PMC like material, with a higher stifihess 
ratio, e.g.. s r = 30) for numerous fiber aspect ratios and two fiber volume fractions, a low (v/ = 0.2) and 
a high (v f = 0.5). The Hashin-Shtrikman-. CMC- and FE-models wore evaluated only at discrete aspect 
ratios: whereas, the VCM-result were computed for many different aspect ratios and thus displayed with a 
solid line. Note, for high stifihess ratios only the lower VCM-bound (which corresponds to the MTM) was 
plotted From Figs. 4 and 5 it is apparent that the influence of fiber aspect ratio on E u is significant (Max. 
error of 35%) for a composite with a high stiffness ratio whereas relatively insignificant (Max. error of 3%) 
for materials with a low stifihess ratio. In both cases the Hashin-Shtrikman estimate lay as expected within 
the upper and lower VCM bounds. The GMC- and FE-results are partially- out of the Mori-Tanka based 
VCM-bounds due to the different inclusion topology (ellipsoidal versus cy lindrical). Generally speaking, 
the GMC-model underestimates the longitudinal stifihess. particularly for high stiflhess ratios, although the 
trend is captured quite nicely . Alternatively, for materials with low stiflhess ratios, good agreement between 
GMC and Hashin-Shtrikman estimates as compared with FE-results are observed given a wide range of 
aspect ratios and fiber volume fractions. 



Figure 4: Comparison of axial modulus ratios for s r = = 3 

The physical reason for the difference between the FE A- results and those of GMC can be inferred from 
Fig. 6, where the longitudinal stress distributions for a load in the a: i -direction are shown. In both FEA 
anal} see shown in Fig. 6, the geometry, boundary conditions, and applied forces are the same. Onl} the 
stifihess ratio between matrix and fiber constituent material is different. Evident!}, these stress contours 
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Figure 5: Comparison of axial modulus ratios for s r = = 30 


indicate that an increase in the stiffness ratio leads to a decrease in the matrix loading and convereelv to 
an increase in the fiber loading (see maximum and minimum stress values). Similarly, the difference of the 
“mean” matrix stress in the volume of material between the two fiber ends as compared with the ‘mean' fiber 
stress becomes greater as the stiffness ratio is increased Considering, the GMC representation, the stresses 
are constant and equal along a longitudinal fiber-matrix row. Therefore the error for lower stiffness ratios 
is decreased, yet when the stiffness ratio increases so too does the error. Figure 6 also gives an explanation 
as too why the GMC stiffness results are slightly under the FE-results in Fig. 4, as the GMC assumptions 
lead to an unloading of the “stiff" fiber and to an additional loading of the “soft" matrix. 

Figures 7 and 8 illustrate the normalized transverse modulus results obtained from the various analytical 
and numerical methods considered Here as one would expect, the fiber aspect ratio has little, if anv, 
effect on the transverse modulus Also, in general the GMC results compare favorably (ma ximum error of 
approximate!} 5%) w ith the FE-results throughout. For both stiffness ratios, the Hashin-Shtrikman estimates 
are equal to the lower VCM bound (which in turn are equal to the original MTM-results). Noteworthy is the 
fact that the GMC-results are much closer to the FE-results for the case of high stiffness ratios, low' aspect 
ratios, and high fiber volume fractions than the other analytical estimates. Similarly, as observed for the 
longitudinal stiffness the qualitative trend of the FE results is better captuied by GMC. 

To help explain this observation, we examine the calculated transverse stress distribution obtained from 
FEA, as shown in Fig. 9. Note, only the case for the low stiffness ratio (s r = 3) is shown in Fig. 9 as 
the FE-results for the higher stiffness ratios appear similar. Relative to the FE-results for the longitudinal 
loading case, the difference in the ‘mean" matrix stress of the volume of material between the two fibers 
in the loading direction and that of the “mean'' fiber stresses for the transverse loading case are small. 
Consequent!}', the inaccuracy caused by having similar stress states in a given row or column of subcells 
within GMC is small; thus explaining the lack of influence the stiffness ratio has on the GMC predictions 
relative to those determined using an FEA unit cell. 

Figures 10 and 1 1 show a comparison of the axial Poisson 's ratios based on the various evaluation methods. 
Again the comparison between FEA and GMC results are quite favorable, irrespective of aspect ratio and 
stiffness ratio; the difference being approximatel} a maximum of 5.2%. Again, the GMC results follow the 
qualitative trend of the FE-results much better than either the Hashin Shtrikman estimates or MTM-results, 
especiallv for the case of high stiffness ratios and low fiber aspect ratios. 
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Figure 6: Example for the longitudinal stresses obtained from FE-ralculation for s r — 3 and .s> = 30 
(V r = 0.2. an = 5.0) 
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Figure 7: Comparison of transverse modulus ratios for s r = 3 
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Figure 8: Comparison of transverse modulus ratios for s r — 30 
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Figure 9: Example of transverse stresses obtained from FE-calculation for .s r = 3 {Vp = 0.2, hr = 5.0) 
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Figure 10: Comparison of axial Poisson ratios for s r = 3 
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Figure 11: Comparison of axial Poisson ratios for s r = 30 

Similarly, results for the in-plane Poisson’s ratio v 23 appear in Figs. 12 and 13. Once again, in general 
the GMC results are slightly above those produced with FEA. The Hashin-Shtrikman estimates being closer 
to the FE-results than the GMC results. The discrepancy between methods being more evident with higher 
fiber volume fraction than higher stiffness ratios. Even so, the maximum error between the GMC predictions 
and those of FEA is at most 4.8%. 

4.2 Elastoplastic Material Behavior 

In this section the elastoplastic stiffness behavior of aligned short-fiber composites will be analyzed using 
primarily the GMC and the finite element method. A comparison of the resulting predictions should confirm 
whether or not the GMC is capable of accurately simulating the elastoplastic stiffness behavior of DRX 
materials. Additionally, results obtained from a semi- analytical elastoplastic shear lag model (see Appendix 
B) are included as well to help illuminate any differences between the methods. The investigated model 
parameters are similar to those of the elastic investigation, with the exception that only two fiber aspect 
ratios (a r — 5 and 25) are examined. Furthermore we restrict ourselves to the case of longitudinal (xj — 
direction) loading so as to allow the inclusion of shear lag results. 

Figure 14 shows the effective longitudinal stress-strain curve produced by the GMC (dotted and dashed 
line), FE unit cell (solid line), and shear lag (symbols) models for the special case of V/ = 0.2, a r — 5, s r — 3. 
Clearly, a large difference in the longitudinal stiffness response exists, depending upon the analysis method 
employed, with the GMC model predictions being significantly softer. The question before us is whether or 
not the lack of shear coupling inherent in the GMC method is the primary cause for the softer elastoplastic 
behavior. To address this question, we introduced a semi-analytical shear-lag-model (and shear-free model) as 
described in Appendix B. A comparison of the shear-lag (open circles) and shear-free (open square symbols) 
models with those of FE and GMC are also shown in Fig. 14. 

Evidently, the 8-subcell (1 fiber subcell surrounded by 7 matrix subcells) GMC results and the shear-free 
model results agree extremely well: whereas, those determined using the elastoplastic shear- lag model agree 
with those coming from the FE unit cell model. Consequently, it may be concluded that: 

1. The difference between the GMC- and FE-simulations is directly related to the lack of normal-shear 
coupling inherent in the GMC formulation. 
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Figure 12: Comparison of transverse Poisson ratios for s r = 3 



Figure 13: Comparison of transverse Poisson ratios for $ r = 30 
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Effective Longitudinal Strain, 

Figure 14: Comparison of FE-, GMC-, Shear-lag-, Shear-free-model results (s r = 3. Vf = 0.2, a r = 5) 

2. The standard (without any modifications) GMC-model Is not well suited for predicting the elastoplastic 
response beharior of DRX materials, as the GMC-model predicts an overly soft plastic response. 

3. The semi-analytical elastoplastic shear-free and shear-lag models introduced in Appendix B appear 
full}* capable of predicting both the GMC- and FEA-results. The primary disadvantage of this shear- 
lag model being its "uniaxial" nature, that is, it can only provide the effective modulus (response) in 
a single direction. 

In addition to the pseudo-square^ fiber (8-suboell) GMC representation (dotted line) in Fig. 14, the 
typical (338 subcell ) circular fiber unit cell (dashed line) is shown in Fig. 14. Clearly, a discretization effect 
is observed and is a result of the piecewise plastification of the matrix subcells between the two fiber ends 
in the case of the finer GMC-model. This effect results in a softer response within the initial plastic range 
(i.e., from 0.2 - 0.5% strain range). After plastification of all matrix subcells the effective hardening slope of 
both GMC unit cells are essentially the same. 

To verify the above conclusions, simulations with other aspect ratios, volume fractions and stiffness ratios 
were conducted with the results being displayed in Figs. 15 - 18. The semi-analytical shear-free model is not 
shown in these figures as its results always agreed with the GMC (pseudo square) simulations. The fitting 
parameter w of the semi -analytical shear-lag-model is found to be only a function of the fiber volume fraction 
i.e. v = v(V/) and takes a value of approximate!} 2 for Vf = 0.2 and approximate!}* 2.5 for Vf = 0.5. The 
dashed lines in all figures represent the effective stress-strain response of the standard circular GMC-model. 
All figures confirm the above conclusions. 

Some characteristic features of the effective finite element stress-strain curves in Figs. 15 - IS are evident; 
for example, increasing either the aspect ratio or volume fraction significantly hardens the composite response, 
particularly for high stiffness ratio materials. The elastoplastic stiffness response may be partitioned into 
three ranges. Figure 19 shows these ranges as a function of effective longitudinal strain. 

Range 1, delineates the domain where linear elastic behavior Is manifested. This is followed by range 2 
in winch a large variation in the hardening slope takes place and range 3 in w T hich approximate!}* a constant 

3 In GMC the geometry of the subcell representing the fiber is square: however in actuality since all fields are taken at the 
centroid of the subcell onh the ratio of subcell to total length is important, no corners are felt. Thus the square fiber is really 
a pseudo-square. 
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Figure 15: Comparison of FE-, GMC- and Shear-lag-model results for s r = 3 and Vf = 0.2 



Figure 16: Comparison of FE-. GMC- and Shear-lag-model results s r = 3 and V/ =0.5 
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Figure 17: Comparison of FE-. CMC- and Shear-lag-model results s r = 30 and Vf- 0.2 



Figure 18: Comparison of FE-, CMC- and Shear-lag-model results s r = 30 and Vf = 0.5 
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Figure 19: Semi-analytical shear-lag model, characteristic ranges and number of yielded truss elements 

(s r = 3, Vf = 0.2, a r = 5) 


effective hardening slope is observed. More detailed investigations into the flow behavior show that at the 
be ginnin g of range 2 the matrix material between the two fiber ends begins to yield first (see Fig. 19, point 
1). After which the whole matrix yields (point 2). During the longitudinal loading step the matrix from the 
side-to-side region tries to move to the end-to-end region. This material flow' leads to a considerable shear 
loading of the side-to-side matrix interphase (between fiber and matrix). In the plastic range large shear 
deformations of the interphase start at the fiber ends (point 3), then grow over the fiber and prevent any 
additional loading of the fiber. In the case of the semi-analytical shear-lag-model this means that one by 
one each non-linear truss element will start to yield. All trusses are plastic at point 4 (see Fig. 19), w hich is 
the beginning of range 3. Note, the stair-stepped dashed line in Fig. 19 indicates the piecewise yielding of, 
altogether, five truss elements. 

Based on these results, the question remaining is whether the GMC can be appropriately modified in 
order to predict accurately the elastoplastic behavior of DRX materials. Considering why the GMC fails, 
it beoomes clear that if the matrix between the two fiber ends yields, the fiber is no longer capable of 
earn ing any additional load beyond the matrix since both the matrix and fiber have the same stress state. 
This continuity of specific stress components along a given row or column of suboells is a direct result of the 
inherent lack of normal and shear coupling with GMC. Given this fact, a simple modification to GMC would 
be the introduction of an artificial matrix material for all subcells within the region between the two fiber 
ends. This new artificial matrix material must be chosen so that the resulting effective stress-strain curve 
produced corresponds with the reference data (be it experimental or analytical). Of course this strategy 
possesses the drawback that the obtained micro-mechanical GMC stress state is no longer based on any 
physical meaning; i.e., fracture, damage or fatigue analyses could no longer be performed confidently w hen 
using the modified GMC-results, as the micro field quantities would be even more questionable than before. 

As a first attempt, the artificial matrix material was assumed to behave as purely elastic, with the 
associated Youngs modulus of the artificial matrix material being that of the matrix itself. Clearly, such an 
assumption results in the overestimation (at least well within the plastic range) of the elastoplastic response 
curve as illustrated in Fig. 20. This response prediction is non-conservative. Ln contrast, it is shown that the 
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Figure 20: Comparison of FE-, GMC-, and GMC-modified simulations (s r — 3, Vf = 0.2, a r =5). 


standard GMC-model provides a lower bound (conservative response prediction) for the effective elastoplastic 
behavior, thus suggesting that an artificial bilinear elastoplastic material might provide a more accurate 
simulation. For such a material definition, four material parameters would be required as input; the}' are 
1) Poisson's ratio, 2) Young s modulus, 3) yield stress, and 4) the plastic hardening slope. With the yield 
stress and hardening slope being the best candidates for modification from that of the matrix material, since 
the elastic range is already accurate!}' predicted, the parameters (E,/y) can be taken to be that of the matrix 
material. Here we restrict ourselves to the modification of the hardening slope alone (keeping the yield 
the same as the of the matrix) and thereby attempt to obtain the unknown artificial hardening slope Eh 
from the consideration of a simple one dimensional spring model (with each spring corresponding to a gi ven 
subcell) which is shown in Fig. 21. 


The following discussion is limited to range 3 in Fig. 19 (i.e., the range where the effective hardening 
slope is approximate!}' constant). Assume an incremental load, AF, is applied to the spring assemblage in 
Fig. 21. This results in an incremental displacement, A u, that depends on both the geometry and material 
properties of the spring assemblage, that is, 


A F = 



_1_ 

+ 


AzLl. 

EZ.A * 


+ =) 

A f ) J 


A u, 


( 2 ) 


where A is the total cross sectional area of the considered unit cell and Af the fiber cross sectional area. 
Em and E* H are the hardening slope and artificial hardening slope, respectively. Alternatively, the global 
longitudinal stiffness behavior can be described by 


AF = 


l 


A u 


( 3 ) 


where E e J^ is the effective longitudinal hardening slope. The combination of Eqs. (2) and (3) , and the 
introduction of the following abbreviations, 
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A F 



Figure 21: Spring assemblage model 


lead to the final equation for the artificial hardening slope: 

E* h S 1.h( 1 - l r) (frM + 7^ ~ *) 

— J * (p) 

H V/(^f-l) + / r (l-<-") 

Note, the only unknown in this equation is the effective longitudinal haidening slope which for example 
may be determined from the semi-analytical shear-lag model (see Appendix B) or obtained from some other 
reference data. 

A comparison of the FE-. GMC-. GMC-bilinear-niodified-model response histories are shown in Figs. 22- 
25. Evidently, as one might expect, the GMC-model, modified by using an artificial bilinear matrix material 
definition in the intermediate subcell, results in a more accurate longitudinal simulation than the previous 
standard GMC-model (see Fig. 14 and that modified using an artificial elastic matrix material definition). 
With this approach, reasonably accurate results can be produced when range 2 Is small, typically for high fiber 
aspect ratios and low stiffness ratios. Alternatively, for situations invoking high stiffness ratio materials, 
low volume fractions and high aspect ratios (e.g., a r = 25, see Fig. 24) simulations using the modified 
GMC approach give the largest error. This Is due to the fact that these conditions induce the smallest area 
(highly loaded) subcell between fibers and thus maximizes the flow between inclusions; thereby producing 
the mast compliant overall response. Clearly, only a small difference exists between FE and GMC results in 
the extensive plastic range. This Is due to the fitting procedure utilized for the simple spring model. Note 
that, again all response histories are conservative when compared w ith FEA simulations. 

More detailed studies were not performed at this time using a more complex artificial matrix material 
definition, since the actual usefulness of such an ad-hoc approach is debatable due to the lack of any real 
physical meaning. Furthermore, a new higher order GMC theory (appropriately termed - high fidelity 
GMC), which provides the necessary normal and shear coupling, has recently been developed by [Aboudi et. 
al.(2001)]. This new formulation, although more computationally demanding than the standard or modified 
GMC (but significant^ faster than the FEA approach) should provide the required balance between speed 
and accuracy when attempting to simulate DRX material behavior. Consequently, further analysis using 
the GMC w ill be reserved until later when the theoretical extension of the present elastic continuously 
reinforced high fidelity GMC version (see [Aboudi et. al.(2001)]) to that of discontinuous and inelastic 
material behavior, Is accomplished 
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Figure 22: Comparison of FE- and GMC-modified model (s r = 3, Vf = 0.2) 



Figure 23: Comparison of FE- and GMC-modified model simulations ( s r — 3, Vf = 
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Figure 24: Comparison of FE- and GMC-modified model simulations (s r = 30, Vf — 0.2) 



Figure 25: Comparison of FE- and GMC-modified model simulations (s r = 30. Vf = 0.5) 
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5 Conclusions 


This paper has focused on evaluating the predictive ability of the generalized method of cells, in the 
context of discontinuous!; reinforced composite (DRX) material behavior; specifically aligned short fiber 
composites. CMC simulations were compared with: i) simulations using other micromechanical based mean 
field models and finite element (FE) unit cell models found in the literature for the elastic range, as well 
as ii) finite element unit cell model and a new semi-analytical elastoplastic shear lag model in the inelastic 
range. Results indicated that the GMC is full} capable of accurately analyzing the elastic behavior of DRX 
materials. However, in the case of inelastic behavior, the standard GMC approach, due to the inherent lack 
of normal-shear coupling, was shown to be unsuitable for predicting the longitudinal elastoplastic behavior 
of aligned short fiber reinforced composites. However, a simple modification of the GMC-model, through 
the introduction of an "artificial" bilinear matrix material definition within the subcell(s) residing between 
inclusion (fiber) ends was shown to improve the overly conservative GMC predictions within the elastoplastic 
regime significantly. This ad-hoc modification to the GMC unit cell definition would allow reasonably 
accurate multiaxial simulations to be made. The case w'hen both the stiffness and fiber aspoct ratios are 
high being the least accurate. Finally, it was shown herein that an elastoplastic shear lag model can 
provide identical uniaxial results to that of the FE unit cell calculations. For truly consistent and accurate 
results it is recommended that a high fidelity micromechanics approach be used when examining inelastic 
discontinuous!} reinforced composite behavior. 
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6 Appendix A: Generalized Method of Cells Formulation 

GMC can predict the elastic and inelastic thermomechanical response of both continuous and discontin- 
uous multiphased composite materials with an arbitrary internal microstructure and reinforcement shape. 
It is a continuum- based micromechanics model that provides closed-form expressions for the macroscopic 
composite response in terms of the properties, size, shape, distribution, and response of the individual con- 
stituents or phases that make up the material'* . These constituent materials can be represented using an}* 
elastic and/or inelastic deformation and life (e.g.. continuum damage mechanics fatigue model) model. The 
periodic nature of composites typically allow s identification of a repeating unit cell that can be used as a 
building block to construct the entire composite. The properties of this unit cell are thus representative of 
the properties of the entire assemblage once the subcell dimensious and appropriate subcell materials are 
selected. Thus, unidirectional long-fiber composites, short-fiber composites, porous materials and laminated 
materials can all be modeled as special cases. 

6.1 Model Description 

Consider a composite material with a periodic structure whose repeating volume element consists of 
N a x N 3 x A/ rectangular parallelepiped subcells. The volume of each one of the subcells is d^hsC ^ , where 
ci, 3 and 7 are running indices: n = 1, . . . , N a ; 3 = 1. . . . N 3 ; = 1. . . . , , in the x\ , xo and x 3 - 

directions, respectively. The total volume of the repeating volume element is dhC where d — X]a=i d 0 * 
h — h >3 1 C = fv In 26, an example is shown for a repeating volume element with A’ a = 2 . 

N 3 = 3 and A r -> = 3. 



Figure 26: Triply periodic repeating unit cell, illustrating subcell dimension nomenclature. 

By an approximate micromechanical anal} sis of the detailed interactions of the subcells of the repeating 
volume element, overall constitutive relations which govern the effective behavior of the multiphase elasto- 
plastic composite, can be established. This analysis relies on the requirements that static equilibrium of the 
materials in the various subcells are ensured, and that continuity of the displacements and tractions between 
neighboring subcells within the repeating volume element, as well as between neighboring repeating volume 
elements is satisfied on an average basis. 

^Xote. each of the subcells can be comprised, in general, by an elastic-viscoplastic temperature -dependent material. 
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As the average behavior of the composite Is sought, it Is sufficient to consider a first order theory in which 
the displacements u^ 3l) in the subcelLs are expanded linearh in terms of the distances from the center of 
each subcell, i.e., in terms of the local coordinates x[ a \ x 2 3 \ and £ 3 ^. Thus, the first order expansion in 
the suboell (a ) is represented as, 

+ x 2 3 > xr 3l) + i - 1, 2, 3 (A-l) 

where are the displacement components at the center of the suboell, and &\ a3 ' l} , X* and 

v •“ 3l ' are micro-variables that characterize the linear dependence of the displacement u i a3 ~" on the local 
coordinates x[ a \ x { 2 3 \ x^K In Eq. (A-l) and the sequel, repeated Greek letters do not imply summation. 
Note that due to the linearity of Eq. (A-l) , static equilibrium of the material within the subcell (a 3~) ) is 
ensured. 

The components of the small strain tensor are given by 


ti 3 ' ] =kdi< a3l) +dju^ ] ) 


ij = 1.2.3 


(A-2) 


where d\ = d/dx* ] , 82 = d/dx [ 2 ] and 8 $ = 8 / 8 x 3 1 . 

The \olume averaged total strains and stresses in the composite are expressed, respectively, as 


Ay A> Ay 


J 




(« .^7) 
ij 


and 


0 = 1 3=1 1=1 


Ay Ns Ay 


«^E E J2 dah3C ^ a ij 


-(>3^) 


(A-3) 


(A-4) 


a=l 3=1 7=1 


with the relationship between the averages of the stress cr\* 3l) , total strain inelastic strain 

and thermal strain e T ^ a3 ^ } = a { - a31 1 AT ( where ' is the coefficient of thermal expansion tensor, and AT 
is the temperature de\iation T — T R from a reference temperature Tr) being given by the classic Hookean 
constitutive equation 


-(o3*>) ^(aSy) /-(a 3y) -I(a^) -T(a3z)\ 


(A-5) 


where ' is the elastic stiffness tensor of the material. 

It has been shown ([Paler and Aboudi (1992)] and [Aboudi(1995)]) that by employing the displacement 
and traction continuity conditions at the interfaces between the subcells of the repeating volume element, 
and at the interfaces between neighboring repeating volume elements, it is passible to eliminate the micro- 
variables and obtain, via a smoothing operation (homogenization), a set of continuum equations that model 
the overall behavior of the multiphase short -fiber composite. This was achieved by establishing relationships 
which connect the microstrains at the subcells to the total overall macrostrains in the composite via the 
appropriate concentration tensors, A and D. 

6.2 Overall Thermomechanical Constitutive Law 

Given these concentration matrices, and D (Q ^\ expressions for the average total strain and 

stress, respectively, in each suboell can be constructed; that is, 


and 


e(“3->) = A {a3l) i + +ef ) 

£.(«*>) = C («.37) r A <a.*)g + D ,a3l) (e{ + ef ) - -I- e T(o37 >)] 


(A- 6 ) 

(A-7) 
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in terms of the uniform o\era!l strain e (the applied macmstrain) and the subcell inelastic and thermal 
strains, where 



and Am involves the elastic properties C la ^ 1 of the subceD material, A c the geometric dimensions of the 
repeating cell only, and J Ls a matrix containing the imposed average (composite) strains. Note, A ici3l) 
and D (0 ^ ) are square matrices of 6 A a A 3 A 7 order. Similarly, notice that the same matrices D io ^' are 
operating on both the inelastic and thermal strains of the subcells. 

It was shown by Dvorak and Benveniste (1992) and Dvorak (1992) that for anv represen t at ive volume 
element under a uniform overall strain c and temperature change AT', which contains a piecewise uniform 
distribution of thermal and inelastic fields (eigenst rains) associated with the applied loading e , the averages 
of the local strains can be expressed in terms of the mechanical and eigenstrain concentration tensors. These 
tensors depend on the local elastic moduli, and on the shape and volume fraction of the phases and are 
therefore constant. Equation (A- 6 ) Ls consistent with this representation of the average subcell strain, w here 
the constant mechanical and eigenstrain concentration tensors are given by A (a5r) and D [a3 *>\ respectively. 

Consequent ly, the following effective elasto plastic thermomechanical law of the composite can be estab- 
lished 

O =B*(e-e , -c T ) (A- 8 ) 

where the effective elastic stif&iess tensor, B*, of the composite Ls given by 


0=1 3=1 7 =1 

and the composite inelastic strain tensor Ls defined as 


i^7) 


(A-9) 


e 1 - 


-B— 

dhC 


A«, A'j AV 




o=l 3=1 7=1 

and the average thermal strain tensor as 


(A-10) 


6 T = 


—B 


<-i 


dhC 


^ 

Y Y Y d <' h3C -' C[a3 ' )} ( Diah )e * - e Tia3 ^ ] )) 

0 = 1 5=1 -1 = 1 


with subcell inelastic, t r s . and thermal, ef . strains being defined as: 

e* = (e /(11,) g/fAiJv.jv,!) 

«r-(« r(i ii),.. J r(A;A ' Ai ) 

The effective coefficient of thermal expansion vector, ct* of the composite, Ls given by 


(A-ll) 


(A-12) 

(A-13) 


a* = e T /AT (A-14) 

The above micromechanically established overaU thermo-inelastic constitutive law (A-S) Ls valid for any 
tv pe of thermomechanical loading (i.e. , any combination of normal, shear and thermal loadings). A significant 
advantage of this constitutive lavv r stems from the fact that it does not rely on any symmetry conditions that 
may exist under certain types of applied loadings. Thus in the implementation of this law, the question as 
to whether such symmetry conditions exist or not is irrelevant. 
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7 Appendix B: Semi-analytical Axisymmetric Elastoplastic Shear 
Lag Model 


Given the elastoplastic investigations under taken in this study, a simple shear-lag model for aligned 
short fiber composites Ls introduced. As mentioned in the background section, shear lag models were the 
first micromechanics models considered for aligned short-fiber composites. Despite some serious theoretical 
flaws (e.g., the ability to predict only the longitudinal modulus E\\ within the elastic range), shear lag models 
have enjoyed enduring popularity; probably due to their algebraic simplicity and their physical appeal. In 
this paper, we have purposefully limited our investigation to the longitudinal elastoplastic stiffness behavior 
so that the proposed modified shear lag model can be appropriate!} utilize. 

Following Cox (1952), the shear lag analyses are focussed on a single fiber of length If and radius ?’/, 
which is encased in a concentric cylindrical shell of radius R. Such a configuration then leads to the key 
assumption of a shear lag model, which Ls that the shear stress between fiber and matrix (r rz , n *p. z. . . 
cylindrical coordinate system) is proportional to the difference in displacement between the fiber surface and 
the outer matrix surface, i.e., 

Tr *^~2jrr} ( w ( R ' z )- w (r f .z)^ (B-l) 

where H is a constant and can be written as 


H = 


2 irG m 

ln(/?/ry)-r 


(B-2) 


It remains to choose the radius R of the outer matrix cy linder. Several choices have been used (see e.g. [Tucker 
and Liang (1999)]). For our purposes it appears most realistic to let R be equal to a radius corresponding 
to the width of the FE or GMC unit cell. The above equations are based on a one-dimensional analysis of 
a fiber surrounded by a matrix layer (i.e., the Poisson effects are neglected (taken to be zero)) and assumes 
linear elastic material behavior. To extend this analysis into the elastoplastic regime of aligned short fiber 
composites, the existing shear-lag model must be modified so as to; 


1. account for elastoplastic material behavior and 

2. the appropriate Poisson effect. 


The elastoplastic material behavior is simply introduced by replacing G m by T the plastic tangent shear 
modulus Gh ( "shear hardening slope") and an incremental formulation of Eq. (B-l), i.e.. 


dr rz (z) = — ^ ' div(R . z) — dw{rf, z )^ 


with 


77 = 2 * Gh 


(B-3) 


(B-4) 


\n{R/r f )-l' 

This new “shear hardening slope" can then be obtained from the P randtl-Pueuss-equ at ion . For a pure 
shear load in the r-z -plane this equation can be written as 


13 dr rz 
de rz = — ■ dr rz + =— r ri 


2G„ 


2 r r *E l 


rz^H 


and leads finally to the shear hardening slope 

dr r ~ 


Gh = 


dt T ~ 


G m E H 

= Vt; — — = const 

Eh + 3G m 


(B-5) 


(B-6) 
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Axial- symmetry 


Matrix 



Figure 27: Sketch of the semi-analj tical shear lag model 


where G h is constant because E H is a constant. Note, within the plastic range a fitting parameter ip Ls 
alar, introduced. This fitting parameter will account for effects which can not be captured by such a simple 
model. 

The Poisson effect is taken into account by considering an axlsymmetric model (shown in Fig. 27 using 
an exploded view) with appropriate boundary conditions on the concentric cy lindrical matrix shell. 

The model posses the following characteristics: 

1. The cross section of the axisvmmetric model is divided into four regions (/, II, III , IV). Further, 
regions III and IV are subdivided into three subregions each, e.g., Ilia , Illb , IIIc and IV a, IVb , 
IV c (see Fig. 27). 

2. Within all regions only normal stress is allowed (i.e., shear stresses are zero) such that rectangular 
cress sections of individual regions remain rectangular after loading. 

3. The longitudinal stress are constant within regions I and II and are piecewise constant within regions 
III and IV (i.e., constant within each subregion). 

4. The points 1-1’, 2-2’, etc. are connected with non-linear truss elements winch characterize the shear 
forces between the fiber and the matrix. 

5. The stiffness behavior of each non-linear truss element, i, Ls based on Eq. (B-l), and follows from 

F; hrusa = T r . 2i TTfli = h(w(R,z) - w(r/,*))z; = C n truss (w(R. z) - w(r f ,z)j (B-7) 

where, e.g., for i = 1, l\ is the distance between points 1 and 2 in Fig. 27 and represents a kind of 
shear-influence length for truss 1. C} russ represents the stiffness of truss 1 between points 1-T. The 
material behavior of the truss element is bilinear and is characterized by a ‘Yield force 1 (which Ls based 
on Tresca’s y ield condition for pure shear within the interface) and can be written as 

K uss.Y = V 2nr f l r 
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as well as by a stiffness hardening slope which can be written as 


Ci russ .H=vHli. (B-9) 

The number of introduced truss elements depends on the desired aspect ratio. For a r = 5 three truss 
elements were used, whereas for a r = 25 six trusses were applied. 

6. Finally, in addition to the symmetry conditions shown in Fig. 27, s}mmetry is also required on the top 
and right outer surface of the axisAmmetric model (see Fig. 27). 

The introduced semi -analytical shear lag model is sohed numerically using the finite element method. 

The so called semi-analytical shear-free model is obtained by assuming that the truss elements are 
characterized by zero stiffness. Such a model represents a standard GMC-model. 
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